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We simulate the homogeneous nucleation of ice from supercooled liquid water at 220 K in the isobaric-isothermal en- 
semble using the MW monatomic water potential. Monte Carlo simulations using umbrella sampling are performed in 
order to determine the nucleation free energy barrier. We find the Gibbs energy profile to be relatively consistent with 
that predicted by classical nucleation theory; the free energy barrier to nucleation was determined to be ~18&b^ and 
the critical nucleus comprised ~85 ice particles. Growth from the supercooled liquid gives clusters that are predom- 
inantly cubic, whilst starting with a pre-formed subcritical nucleus of cubic or hexagonal ice results in the growth of 
predominantly that phase of ice only. 

PACS numbers: 64.60.Q-, 64.70.D-, 82.60.Nh, 64.60.qe 



I. INTRODUCTION 

It is not unreasonable to expect water to freeze into ice 
when it is cooled below its freezing point of 0°C in normal 
conditions. However, despite everyday experience, it is pos- 
sible to cool very pure water droplets down to approximately 
-45 °C at standard pressure without observing freezing. 1-3 In 
the absence of an external nucleation seed, water can remain 
in its supercooled state for exceedingly long periods of time. 
Classical nucleation theory (CNT) predicts the existence of 
a kinetic barrier to nucleation due to a competition between 
the unfavourable interfacial free energy and the favourable 
bulk free energy of growing clusters, rendering homogeneous 
nucleation extremely slow, as a cluster of some critical size 
must spontaneously form before large-scale crystallisation can 
occur. ~ In other words, homogeneous nucleation is a rare 
event, where the average time between nucleation events is 
many orders of magnitude greater than the duration of the 
event itself. 

Although in practice, heterogeneous nucleation of ice dom- 
inates the freezing process, homogeneous nucleation is inter- 
esting not only as a starting point towards more practical re- 
sults, but also because it is thought to be important in pro- 
cesses such as the formation of ice in clouds, governing both 
global climate change and short-term weather, 7-1 1 and the cry- 
opreservation of biological tissues. 12 On the other hand, the 
fact that homogeneous nucleation is a very rare event (or, 
equivalently, that supercooled water is metastable with respect 
to ice in the absence of external nucleants) is likewise crucial 
in some phenomena, where it is important that water remain 
a (supercooled) liquid, such as electrification and rainfall in 
clouds. 13-16 

As nucleation is a rare event, it is intrinsically difficult to 
study. There are two principal factors potentially govern- 
ing the slowness of nucleation: the nucleation free energy 
barrier could be very large or the dynamics of the process, 
excluding the kinetic barrier to nucleation, might be inher- 
ently slow. Indeed, the highly ordered hydrogen-bonded liq- 
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uid structure, especially at lower temperatures, slows diffusion 
considerably; so much so, in fact, that it has been suggested 
that the Gibbs energy of diffusion activation corresponding to 
the transfer of water molecules across the ice-water interface 
may be rate-determining. 2 ' 14 ' 17 However, other experimental 
work suggests that the temperature dependence of the nucle- 
ation rate primarily reflects the behaviour of the free energy 
barrier of nucleation rather than diffusion activation. 18 

At atmospheric pressure below the freezing point, the ther- 
modynamic form of water is hexagonal ice (ice Ih) with the 
lonsdaleite structure; an important metastable phase in these 
conditions is cubic ice (ice I c ) with the diamond structure in 
the oxygen atom arrangement. Both are proton disordered and 
obey the Bernal-Fowler ice rules. 19-21 Cubic ice appears to 
form as the predominant (if metastable) phase between 130K 
and 150 K, but is converted into hexagonal ice above ap- 
proximately 200 K. It is thought to play an important role in 
certain atmospheric processes. 22 ^' 5 Although cubic ice is al- 
ways the metastable form, the energy released on conversion 
into hexagonal ice is only between 13 and 50Jmol -1 . 21 By 
contrast, for the TIP4P model of water, 26 cubic ice may, de- 
pending on the precise conditions of the simulation, be more 
stable than hexagonal ice, although the energy difference 
may be less than the measurement uncertainty. Both forms 
of ice have been observed in nucleation experiments in various 
conditions. 18 ' 22 ' 29 ' 30 

Ice crystallisation has been studied in a significant number 
of simulations using all-atom models of water; ' 1-47 however, 
despite the obvious importance of the process of ice nucle- 
ation, there have not been many successful attempts at simu- 
lating the process so far. Crystallisation has been observed in 
an individual molecular dynamics (MD) trajectory at a fixed 
density lower than the liquid state density, 34 and it remains the 
only successful brute-force simulation of ice nucleation with 
an all-atom model. Nucleation under 'special' conditions has 
been observed, such as in simulations of very small numbers 
of particles in a strong electric field 31 ' 32 or near a surface; 39 ' 44 
however, in these studies, the ice clusters that grow quickly 
span the simulation box. An alternative approach is to use 
rare event methods to compute the free energy landscapes for 
crystallisation. 
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Such free energy landscapes have recently been calcu- 
lated in several studies. Radhakrishnan and Trout used two- 
dimensional umbrella sampling, 36,37 Quigley and Rodger used 
a metadynamics approach, and Brukhno and co-workers 
used replica exchange umbrella sampling. The first two stud- 
ies mentioned use 'global' rather than 'local' order parame- 
ters. Although these order parameters detect the increase in 
order as crystallisation occurs, they do not specifically refer 
to the formation of a local crystal nucleus, and so it is not to- 
tally clear to what the variables in the energy landscape phys- 
ically refer. Furthermore, there are grounds to think that the 
pathways that these order parameters enhance may not be rep- 
resentative of real ice nucleation. An increase in these order 
parameters is dependent on the accumulation of local Stein- 
hardt vectors such that they do not cancel out. Once a suffi- 
cient random fluctuation occurs in one of the orientations, it is 
driven onwards, which ensures that all the ice-like molecules 
are orientated in a correlated fashion; indeed, the orienta- 
tional coherence affects even those water molecules that re- 
main liquid-like. This appears to be at odds with a local pic- 
ture of nucleation, since a molecule's behaviour is not just in- 
fluenced by its local environment, but also, in a non-physical 
way, by the orientation of (ice) molecules distant from it. The 
use of global order parameters to drive crystallisation has also 
been shown to induce entropic break-up of small clusters. 49 

The study by Brukhno and co-workers used novel local 
order parameters based on a maximum projection director 
approach, which do allow the calculation of a free energy 
profile as a function of an order parameter measuring the in- 
creasing crystallinity of a nucleus. However, because their or- 
der parameter is not rotationally invariant, but induces crys- 
tallisation with a specific orientation with respect to the sim- 
ulation box, this bias is likely to lead to a significant error in 
the free energy barrier, because it excludes many other pos- 
sible pathways (which could have a lower free energy). The 
space-fixed nature of the order parameter also induces a non- 
local and non-physical orientational coherence in the growing 
ice cluster. 

The ice form generated in these simulations has varied from 
pure cubic ice, to pure hexagonal ice, and to mixtures with 
predominantly hexagonal ice. The free energy barriers com- 
puted in these studies are intriguing. All these simulations 
used the TIP4P water model and presented results at a temper- 
ature of 1 80 K and at standard pressure, which represents an 
approximately 22% supercooling. Brukhno and co-workers 
find a free energy barrier of approximately 130 k^T, although 
this value should be largely discounted because of their use 
of a non-rotationally-invariant order parameter, and, as they 
themselves note, because their system was not fully equili- 
brated. 

The barriers reported by Radhakrishnan and Trout (63 &b T) 
and Quigley and Rodger (19k^T) are more interesting. If 
the equilibrium free energy landscapes have been obtained, 
these values should be lower bounds to the free energy bar- 
rier one would get for the 'perfect' order parameter because of 
configurations at the top of the barrier that are not truly rep- 
resentative of the intermediates between the two states. As 
one would generally expect an order parameter that follows 
the largest crystalline cluster within the system to provide a 
better representation of the intermediates, one would expect 



the free energy barrier obtained with such an approach to be 
larger than that obtained using global orientational order pa- 
rameters. Classical nucleation theory can provide an estimate 
of the former; the CNT expression for the Gibbs energy bar- 
rier of a spherical cluster comprising N particles is AG(N) — 
—NAfusH + y[36n(N/p) 2 ] 1 ^ 3 , where p is the number density 
of the crystalline phase, Af USJ u is the change in chemical po- 
tential on fusion and y is the interfacial free energy between 
the two phases. An estimation of the classical nucleation the- 
ory result for TIP4P water at 180K gives a free energy barrier 
of AG « 35&b^ and a critical nucleus size of 106. To ob- 
tain these, we use the simple approximation that Af us /i(T) « 
Af us #m((7fus — T)/Tf us ), and we insert the appropriate values 
for the TIP4P model. We note that one generally expects the 
interfacial free energy to decrease with decreasing tempera- 
ture, as has been inferred from experiment, ' and so the use 
of y from coexistence likely leads to an overestimation of this 
barrier. 

One is tempted to ask why there is such a big difference in 
the barrier, and in the opposite direction to what is expected. 
The first option is that classical nucleation theory's shortcom- 
ings are responsible for a severe underestimation of the free 
energy barrier; however, although it is well-known that there 
are deficiencies in the framework of CNT, such a big differ- 
ence nevertheless seems surprising. Another option might be 
that the free energy landscapes computed in simulations are 
not fully at equilibrium; to find the free energy barrier of nu- 
cleation, not only the end points, but the entire path space must 
be equilibrated. The latter may have been a problem in the pre- 
vious simulations, as the nature of the order parameters may 
bias the system to locate more easily pathways that are high 
in free energy, but where the dynamics along the reaction co- 
ordinate (once the free energy barrier has been negated) are 
more computationally tractable to observe. For example, the 
coherence induced by the use of global order parameters may 
make it easier for water molecules to join the ice crystal if they 
can 'feel' the orientation of the crystal, even if this is not the 
natural mode of growth. 

Given the importance of water, it is perhaps surprising that 
free energy profiles for ice nucleation using a local, rotation- 
ally invariant order parameter have not been computed for an 
all-atom model of water. The problem seems to be that the nat- 
ural dynamics associated with ice crystallisation are extremely 
slow. For example, even the growth by a few layers of a planar 
ice- water interface below the freezing point, which is a purely 
downhill process in free energy, can take weeks of computer 
time with a TIP4P-type model. 38 The fastest rate of growth 
of such an interface is just below the freezing point, 45-47 but 
this is where the free energy barrier to nucleation is very large. 
Intriguingly, it has been shown that ice crystal growth occurs 
more rapidly in molecular dynamics simulations than in corre- 
sponding Monte Carlo ones, which is not typically the case 
in simulations of other materials. One possible rationalisation 
for this is that collective motion possible in molecular dynam- 
ics simulations helps to speed up the dynamics of cluster reor- 
ganisation in the process of crystallisation. 

Here, we present simulations with a somewhat more lim- 
ited aim, but which we hope present a step towards under- 
standing the homogeneous nucleation of ice. We wish to look 
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at a model of water that does not have the full complexity 
of TIP4P-type models, namely the monatomic model of wa- 
ter (MW) proposed by Moore and Molinero, 57 which has al- 
ready been used to study both the supercooled liquid 58 ' 59 and 
ice crystallisation (in the bulk 60-62 and in confinement 63 ' 64 ). 
Furthermore, some nucleation rates and critical cluster sizes 
have very recently been reported from molecular dynam- 
ics simulations near the liquid transformation temperature of 
202 K and in the range 220 K to 240 K. This potential is 
essentially a modification of the three-body Stillinger- Weber 
potential for silicon with a greater weight given to tetrahe- 
drality. Despite being monatomic, it provides a surprisingly 
good structural and thermodynamic representation of water; 
however, unlike all-atom water models, it crystallises rela- 
tively readily. 60 The absence of explicit hydrogens appears to 
change dramatically the rate of nucleation. This could be be- 
cause of a change in the free energy barrier to nucleation, or 
because the dynamics of the process are much quicker. In- 
deed, it is the unnaturally fast diffusion in the model which is 
most at odds with experiment. A possible rationalisation of 
this behaviour may be that the complication of hydrogen bond 
flipping is no longer an issue. For example, if one molecule 
were to flip in TIP4P water, to continue obeying the ice rules, 
a whole series of water molecules would need to flip, which 
would clearly be an unlikely event. When hydrogens are re- 
moved, this complication no longer exists. In this work, we 
make use of the model's enhanced crystallisability to compute 
the free energy profile for nucleation in the MW model and 
hence to understand better the thermodynamics of ice nucle- 
ation. 



II. SIMULATION METHODS 

To simulate the system, we use the standard Metropo- 
lis Monte Carlo approach. In addition, we use the um- 
brella sampling ' method with windowing in the isobaric- 
isothermal ensemble to calculate the free energy profile. In 
the umbrella sampling simulations, we iteratively update the 
weights depending on the number of steps spent at each value 
of the order parameter in the previous iteration. To test if the 
simulation is equilibrated, we need to ensure (a) that there is 
(approximately) even sampling throughout the window, and 
(b) that there is sufficient exchange between order parameter 
values. A further check to ensure that we have sufficiently 
equilibrated each window is to ascertain that the potential en- 
thalpy as a function of the order parameter is consistent in the 
overlapping regions. 

In simulations, care must be taken to ensure that the ice 
clusters do not begin to span the box. Such configurations are 
a non-physical consequence of the periodic boundary condi- 
tions and should be excluded from our sampling of the nucle- 
ation free energy profile, but can arise due to fluctuations to 
a more cylindrical shape. To reduce the probability of span- 
ning clusters arising, we simulate the nucleation process in a 
cubic box, as if we allow box dimensions to vary indepen- 
dently, random fluctuations in the liquid phase often produce 
long cuboid boxes in which the formation of spanning clusters 
is more likely and which can thus crystallise readily, albeit in 
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FIG. 1. A typical probability distribution for all pairs of ^(z, j) = 
93(0 " Q?>U)i where particles i and j are within 3.6 A of each other. 
The three states depicted have all been equilibrated at 220 K and the 
ice structures are not, therefore, 'perfect' . This figure is analogous to 
those in Refs and ; we reproduce it here for convenience. 



a non-physical manner. We should note that the spanning of 
clusters only happens very occasionally in our simulations be- 
cause we use a sufficiently large number of particles that for 
the size of the crystalline clusters we consider, their surface to 
volume ratio is minimised by a spherical rather than a cylin- 
drical shape. 

The order parameter we use is the size of the largest ice 
cluster in the system. The choice of an ice classification pa- 
rameter is not trivial, as it must on the one hand be local 
and be able to induce the growth of a small ice cluster, but 
on the other hand, it must be strict enough to encourage the 
growth of ice, rather than a structure only locally reminiscent 
of ice, but with no long-range order. We use Steinhardt- style 
local parameters 48 ' 49 similar to those used previously in stud- 
ies of tetrahedral liquids; 64 ' 68 namely, we use the (2/ + 1)- 
dimensional real vector order parameter <//(/) with compo- 
nents m G [—/,/] (m G Z) given by 



qim(i) 



1 



^neighs (0 



AWghs (0 j 



S lm(0ij, <Pij) , 



(1) 



where S/ m are the spherical harmonics, Af ne i g h s is the number 
of neighbours of particle i, 6 and (p are the polar angles, and 
we use 1 = 3. For an up to five-fold decrease in computation 
time, we use real spherical harmonics 69 as opposed to complex 
ones. The sum in j is over all neighbours within 3.6 A. We 
then calculate j) = qi(i) ■ qi(j). The resulting quantity 
ranges between —1 and +1, with cubic ice exhibiting a peak 
at —1 and hexagonal ice exhibiting a large peak at —1 (cor- 
responding to staggered bonds) and a small peak near —0.1 
(corresponding to eclipsed bonds). By plotting the distribu- 
tion of J3 values for the liquid and the solid phases (Fig. 1), 
a critical threshold can be determined as the first point where 
the probability of being in the solid phase is zero. This scheme 
works very well with fully-formed crystals and liquid water, 
which are thus excellently differentiated; however, if we used 
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this as our classification parameter in attempting to grow ice 
clusters, this would unfairly bias us towards cubic ice, as we 
would accept imperfectly formed cubic ice (with 3 connec- 
tions) as ice-like, but imperfectly formed hexagonal ice with 
three suitable connections would not be accepted as ice-like if 
two of these connections were in the < —0.82 region and 
one were in the vicinity of —0.1. If we wish to ensure that 
we do not bias simulations towards cubic ice, we need to be 
careful about our classification definitions. Our overall clas- 
sification parameter is thus defined as 

^neighs (0 

"con(0= £ r{d 3 (ij)), (2) 
7=1 



where 




1 if [(* < -0.82) A (-0.145 < x < -0.065)] 



(3) 

We classify a particle as ice-like if n con > 3 and as liquid- 
like otherwise. This gives perfect identification in both equi- 
librated cubic and equilibrated hexagonal ice. Although it ap- 
pears that liquid water, having a considerable proportion of 
particles in the region where ice Ih has its second peak, 
might be mistakenly classified as ice, it is in fact almost 
nowhere that it has a sufficient number of connections with 
this value of J3, and we find a misidentification rate in super- 
cooled liquid water of approximately 0.8%, suggesting that 
this is an excellent classification parameter. The order param- 
eter we use in umbrella sampling simulations to track nucle- 
ation is the size of the largest cluster of particles classified as 
ice, where two particles belong to the same crystalline cluster 
if they are both ice-like and within 3.6 A of each other. 

To calculate nucleation rates, we perform MD simula- 
tions using the LAMMPS molecular dynamics code in a 
Bennett-Chandler- type approach. 5,7 1-73 We assume that the 
crystal nucleation rate R is related to the Gibbs energy bar- 
rier AG through an Arrhenius-type equation R = KP(n C n t ), 
where n cr i t is the size of the critical cluster, P(n C n t ) = 
exp(— AG(n cr i t ) / faT) and K is the kinetic pre-factor. We 
can assume 71 that the kinetic pre-factor can be written as 
k = Zp Xiq f ncdi , where Z 2 = |AG"(rc crit )| /2nk B T is the Zel- 
dovich factor, f Hcrit is the rate at which particles are attached to 
the critical cluster, and pn q is the number density of the super- 
cooled liquid. The attachment rate f HcTit can be expressed as 
the diffusion of the cluster size at the top of the barrier by a 
modified Einstein relation, /„ crit = lim^oo ((n(t) — ft C rit) 2 ) /2f. 
Although this is considerably simpler than the full Bennett- 
Chandler scheme, 74,75 it has been shown to agree very well 
with other methods of calculating rates in crystal nucleation, 
such as forward flux sampling and direct measurement from 
molecular dynamics. 73 

III. RESULTS 

A. Free energy profile 

In our Monte Carlo simulations, we simulated 1400 MW 
particles at lbar pressure and 220 K, which corresponds to 




Size of largest cluster 



FIG. 2. The free energy profile of MW nucleation as a function of the 
size of the largest crystalline cluster in the system. Simulation results 
from different windows are depicted in alternating colours to show 
their overlap. The dashed line corresponds to the classical nucleation 
theory prediction; dotted lines depict fits to the simulation data. The 
free energy profiles for ice nucleation seeded with hexagonal crystal 
clusters (Ih) and ice nucleation directly from the supercooled liquid 
(FG) are shown. 



a 20% supercooling. The umbrella sampling used partially 
overlapping windows of various sizes until results were con- 
sistent and the criteria for equilibration given in Section II 
were deemed to have been fulfilled. On average, this in- 
volved approximately 5 x 10 10 Monte Carlo steps for each 
window. The simulations involved three distinct scenarios: in 
one case, crystal clusters were grown directly from the super- 
cooled liquid; in the other two cases, we started simulations in 
the initial few umbrella sampling windows with seed clusters 
of suitably equilibrated hexagonal and cubic ice, respectively. 
These were then allowed to shrink and grow, and clusters that 
grew into overlapping window regions were taken as initial 
starting points for the umbrella sampling simulations in later 
windows. We must emphasise that the cubic and hexagonal 
ice simulations were not constrained in any way to form a 
particular phase of ice; we merely introduced small clusters 
taken from a particular crystal structure into the liquid, and 
the growth and shrinkage proceeded unconstrained from those 
starting configurations. 76 The free energy profiles correspond- 
ing to crystal nucleus growth from the supercooled liquid and 
from equilibrated hexagonal clusters are shown in Fig. 2; for 
clarity, we have not plotted the results for nucleation starting 
from cubic clusters, as the results are essentially identical to 
the hexagonal ice case. 

Some snapshots of the nucleation process are shown in 
Fig. 3. Interestingly, the type of ice that grows during this 
nucleation process depends on the starting point. If we be- 
gin a simulation from the supercooled liquid, the ice clusters 
that form are initially mixed phases of hexagonal and cubic 
ice. As these 'mixed' clusters grow, however, they become 
predominantly cubic (except for surface particles, for which a 
classification as cubic or hexagonal is less well-defined 77 ); the 
proportion of core ice particles classified as being hexag- 
onal drops from almost unity to approximately 10% by the 
time the crystal nucleus comprises 50 particles (Fig. 4). Start- 
ing from pre-formed cubic clusters in the lower umbrella sam- 



FIG. 3. Representative nucleation snapshots from umbrella sampling simulations. Particles within 3.6 A are connected with lines; the lines 
are yellow within the largest ice cluster and cyan elsewhere (in the top row only). Particles are not connected across the periodic boundary. 
Spheres represent particles classified as ice: red spheres correspond to cubic ice, orange spheres correspond to hexagonal ice and pink spheres 
(in the top row only) correspond to ice particles not within the largest crystalline cluster. In (a), a 30-particle cluster as nucleated from the 
supercooled liquid is shown. In (b), an 83-particle cluster as grown in a simulation initially seeded with an equilibrated cubic ice cluster is 
shown. In (c), a 165-particle cluster as grown in a simulation initially seeded with an equilibrated hexagonal ice cluster is shown. In each case, 
the top and bottom pictures depict the same cluster from different perspectives; one within the liquid framework and one showing solely the 
largest crystalline cluster. 



pling windows also results in primarily cubic ice growth. By 
contrast, starting the simulation with a pre-formed hexagonal 
cluster close in size to, but slightly smaller than, the critical 
size, results in essentially pure hexagonal ice growth. This 
suggests that in this model, ice prefers to grow in the phase of 
the underlying crystal nucleus. We have not been able to find 
a significant difference in the nucleation free energy profiles 
corresponding to cubic and hexagonal ice. However, our simu- 
lations show that the free energy barrier to nucleation of freely 
grown (mainly cubic) ice is ~2Ak^T compared to ~18&b^ 
for pre-formed (hexagonal or cubic) ice, and the critical size is 
^1 14 compared to ^85 (see Fig. 2). The primary difference 
arises in the early stages of nucleation, where the ice grown 
directly from the supercooled liquid relaxes from mixed clus- 
ters to a much purer form of cubic ice. Nevertheless, the 'cu- 
bic' ice grown from a supercooled liquid has more hexagonal 
defects than the corresponding cubic ice grown from a small 
perfectly cubic crystal nucleus. 

In a very recent study 62 of homogeneous ice nucleation in 
the MW model, Li and co-workers studied the kinetics of nu- 
cleation using forward flux sampling at a range of tempera- 
tures, including at 220 K. They found a critical crystal nucleus 
size of 265, which seems at first glance to be a considerably 
larger cluster than the sizes we obtain. However, the order pa- 
rameter used in their study was different from the one we use; 
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FIG. 4. The proportion of core ice particles classified as hexagonal 
for the set of simulations in which the crystalline cluster was grown 
directly from the supercooled liquid. Error bars show the standard 
deviation for the population of configurations at each ice cluster size. 



the principal difference is that all the neighbours of their clas- 
sified ice particles were counted as being ice-like in their work 
(called 'surface' particles), whilst we treat such particles as 
liquid-like. We calculated the sizes of our critical clusters us- 
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ing their order parameters (utilising / = 6 spherical harmonics, 
a real-space neighbour cutoff of 3.2 A and a dot product cutoff 
of 0.5). When surface particles are not included, the cluster 
size is approximately equal to the values we obtained; how- 
ever, including surface particles gives critical sizes of approx- 
imately 150 and 200 (corresponding to sizes of 85 and 114 
without surface particles, respectively). Although this is still 
less than their reported 265, the disagreement is not nearly as 
significant as might initially appear. We can rationalise the dis- 
crepancy by noting that our system is able to relax locally due 
to the equilibration afforded it by umbrella sampling, whilst 
their system is driven forward aggressively by the forward flux 
sampling algorithm. 

Similarly, Moore and Molinero have recently performed 
crystallisation simulations near the liquid transformation 
temperature; they estimate that the critical nucleus at 208 K 
comprises 120 particles, whereas our simulations suggest, fol- 
lowing histogram reweighting, that the critical nucleus at that 
temperature contains approximately 50 particles for the sim- 
ulations grown directly from the supercooled liquid. As with 
the study by Li and co-workers, the discrepancy arises pri- 
marily from Moore and Molinero 's inclusion of what they call 
'intermediate ice' in their largest cluster. 

The question of whether 'surface' particles are ice-like or 
liquid-like is not easy to answer, and whether restrictive or 
generous order parameters are more realistic is unclear; this 
is indeed a rather general problem affecting all nucleation 
studies. 79 It is certainly the case in our experience that parti- 
cles near the surface of 'core' ice can be very neatly arranged 
and thus reminiscent of the ice phase; but equally, there are 
many examples of particles adjacent to the ice cluster that have 
a structure totally incompatible with ice (such as having more 
than five neighbours). Nevertheless, although the choice of 
order parameter will affect the appearance of the free energy 
landscape (and the sampling efficiency), it should in princi- 
ple have little effect on the barrier height and no effect on the 
rate of nucleation, provided that it is still a good parameter to 
describe the nucleation process. 

Moore and Molinero studied the crystallisation of super- 
cooled MW model water at 180K, where they observed ice 
formed from the supercooled liquid to consist of intercalated 
layers of cubic and hexagonal ice in the approximate ratio 
2:1. ' We have not seen such behaviour in our simulations; 
an obvious possible reason for this difference is that their sim- 
ulations were run at a lower temperature: indeed, they were 
run below the liquid transformation temperature of 202 K, 
where the structure of the liquid phase changes. However, Li 
and co-workers observed the ice grown in their simulations to 
be initially hexagonal, but transforming to a mixture of cubic 
and hexagonal ice in the approximate ratio 1 : 1 in later stages 
of the nucleation at temperatures including 220 K. Umbrella 
sampling allows clusters both to grow and to shrink and thus 
to equilibrate locally, whereas their observations may be a re- 
sult of kinetics dominating the growth pathway. In this regard, 
it is noteworthy that we have occasionally observed a hexago- 
nal defect on cubic ice leading to further hexagonal ice growth 
and vice versa, and faster growth encourages more defects. 

It is also interesting to note that the system can become 
trapped in different regions of path space. This is the case 
not only for 'pure' cubic and hexagonal ice clusters, for which 



we would expect it to be difficult to interconvert once growth 
has started, but also between the ice clusters grown directly 
from the supercooled liquid, which were predominantly cubic, 
and the lower free energy pure cubic ice clusters. Such non- 
equilibrium effects are not new in studies of nucleation; for 
example, analogous behaviour has been observed in the crys- 
tal nucleation of a binary suspension of oppositely charged 
colloids. 81 

From classical nucleation theory, we can estimate the free 
energy of nucleation using the formula given in the Introduc- 
tion. We can estimate Af usj u by performing Gibbs-Helmholtz 
integration from the coexistence point (T\ — ). To do 

this, we can evaluate 

= constant -f^+f^dr, (4) 

where (U) m is the mean molar potential energy and (V) m is 
the mean molar volume measured in a simulation, averaged 
over a reasonably large number of Monte Carlo cycles. We 
perform such simulations over a range of temperatures for 
both the liquid state and the hexagonal ice state and then find 
best-fit equations for the mean potential energy and the mean 
volume, and integrate appropriately. Using this approach in a 
series of short simulations, we find that the chemical poten- 
tial difference between the two phases is Af us /i(220K)/&B ~ 
1 1 8 K per particle, in surprisingly good agreement, consider- 
ing the degree of supercooling, with the value obtained from 
the simple approximation that Af usj u(r) « Af us if m ((7f us — 
T)/Tf us ), 7 where we use Af us H = 1.26kcalmol _1 and 7f us = 
274.6K, 57 giving A fus ^(220K)/^ B ~ 126 K per particle. We 
are not aware of any calculated values of the interfacial free 
energy term for the MW potential and so we use that de- 
termined for the TIP4P potential, 24mJm~ 2 . 54 The classi- 
cal nucleation theory free energy profile is shown for com- 
parison in Fig. 2; it is qualitatively similar to the computed 
profiles, although the barrier is somewhat higher. One as- 
pect to bear in mind is that the CNT profile passes through 
the origin, whereas for our order parameter, the size of the 
largest crystalline cluster has an average value of two in the 
supercooled water. Interestingly, a line of best fit corre- 
sponding to AG(N) = -A ius nN + t 2 N 2 / 3 +t 3 N 1 / 3 +f 4 , where 
ti are parameters determined by regression (AG(N)/k^T = 
-0.547V + 3.457V 2 / 3 + 0.847V 1 / 3 -6.13), suggests the interfa- 
cial free energy attributable to the lowest free energy pathway 
is 7~ 22.9mJm -2 (also shown in Fig. 2). 82 A similar analy- 
sis for the system grown directly from the supercooled liquid, 
where we fit all four parameters, including Af us /i, suggests 
that A fus ^(220K)/^ B ~ 120K and y w 26.2mJm~ 2 ; in other 
words, the chemical potential does not change significantly, 
but the interfacial energy becomes somewhat less favourable 
for 'mixed' crystal growth. Given the physically reasonable 
values obtained from these fits, we can conclude that classical 
nucleation theory appears to apply to a reasonable degree to 
the nucleation of MW ice; an accurate numerical determina- 
tion of the interfacial free energy for the MW model would 
help in assessing this observation. It is also interesting to note 
that the enthalpy associated with the growth of the crystal nu- 
cleus is a monotonically downhill function of the nucleus size; 
this result implies that the barrier to nucleation is primarily en- 
tropic in nature. 
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Li and co-workers estimated the interfacial free energy from 
a fit of the nucleation rate to the CNT result as a function of 
the temperature, assuming that 7 is temperature-independent 
and that the critical nucleus is spherical. They reported a 
value of 7 = 31.01 mJm -2 , which is higher than what we 
have determined. In order to compare our results with theirs, 
we have calculated the order parameter used by Li and co- 
workers for a large number (~60 000) of independent config- 
urations taken along the nucleation pathway for the simula- 
tions started from the supercooled liquid; we weighted each 
value with the associated free energy from Fig. 2. This pro- 
duced an approximate free energy profile corresponding to nu- 
cleation measured with their order parameter. Fitting a CNT- 
like expression, AG(N)/k B T = t { N + t 2 N 2 / 3 +t 3 N 1 / 3 +t 4 , to 
this free energy profile allows us to calculate the interfacial 
free energy associated with the process. If we fix t\ to the bulk 
chemical potential change (t\ = — Af USJ u(220K)/^B^) ? we find 
that 7~ 34.9 mJm -2 , in reasonable agreement with the result 
of Li and co-workers. However, this fit gives t 3 = —8.9 and 
t\ — — 0.04; this is a much larger value of t 3 than we obtained 
with our order parameter, and because t-$ does not appear in 
the CNT expression, this suggests that this may not be the 
best fit to the data. By contrast, if we fit all four t{ parameters 
and interpret them using classical nucleation theory, we find 
that A fusJ u(220K)//;B w 81.52K and 7^ 20.8mJm- 2 ; here, 
t3 = 0.28, suggesting that the fit to classical nucleation theory 
is much better. 

We suggest that this fitting process may provide a measure 
of how reasonable the order parameter used to track the nucle- 
ation process is; an unconstrained fit should reproduce Af usj u 
reasonably when fitting data to CNT. This is the case when us- 
ing our order parameter, but not when using the order parame- 
ter of Li and co-workers; their order parameter is too generous 
and the result is a chemical potential difference that does not 
appear to favour ice as much as it should, and a consequently 
lower interfacial free energy, since much of the ice cluster is 
liquid-like. If the assumption is that Af us /i can be approxi- 
mated by its bulk value in this case, then the interfacial free 
energy will necessarily be overestimated in compensation. We 
suggest, therefore, that counting 'surface' particles as being 
part of the ice cluster may result in an order parameter that is 
not well suited to being utilised in a fit to classical nucleation 
theory. 

Since the MW potential is a reparameterisation of the 
Stillinger- Weber potential for silicon, it is worthwhile also 
briefly to compare the free energy profile to that calculated 
for the nucleation of Stillinger- Weber silicon, which Beaucage 
and Mousseau simulated in the canonical ensemble at 25 % su- 
percooling and for which they determined the critical nucleus 
size to be approximately 175. Although the conditions and 
the order parameter used to track the nucleation process dif- 
fer from our simulations, we note that the increased weight 
given to tetrahedrality in the MW potential appears to result in 
a smaller critical nucleus size. This is not unexpected in the 
light of classical nucleation theory, as the greater similarity of 
the liquid to the solid presumably decreases the interfacial free 
energy. 

From a visual inspection of configurations, the crystalline 
clusters seemed to be fairly spherical. To quantify this ob- 
servation, we have computed two measures of the sphericity 
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FIG. 5. A plot of two sphericity parameters against the size of the 
largest crystalline cluster calculated for snapshots of the system along 
the local order parameter used to drive nucleation in our system. Er- 
ror bars show the standard deviation for the population of configura- 
tions at each cluster size. The diagram shows results for nucleation 
from a supercooled liquid; the other two systems behave analogously. 

of the clusters. The first is based upon the moment of inertia 
tensor, which has elements 

N 

hj = m k( r k 5 ij ~ r kirkj), (5) 
k=\ 

where is the z-th component of the vector between the cen- 
tre of mass and particle k, is the overall magnitude of this 
vector and fyj is the Kronecker delta; i and j refer to cartesian 
co-ordinates in an arbitrary laboratory frame of reference. The 
eigenvalues of this tensor give the principal components of the 
moment of inertia; we use these to quantify the orientational 
sphericity using the asphericity parameter 

_ (jg - Iyy) 2 + (jg ~ kz) 2 + (Iyy ~ kz) 2 _ 

2(I XX +Iyy+I ZZ ) 2 ' ^ 

which ranges between zero for a perfectly spherical cluster 
and unity for extremely elongated ones. Any cluster with s ;§ 
0. 1 appears very much spherical by visual inspection. As the 
clusters grow in our simulations, their sphericity increases, as 
can be seen in Fig. 5, although this is to some extent a result 
of the discrete nature of clusters being felt more at smaller 
sizes. We also characterise sphericity in terms of the radius of 
gyration, namely using 58 

* g = 1^73 " h where = h %%} ri ~ Tj ' 2 ' (7) 

c is the size of the cluster and the 1.5c 1 / 3 factor reflects the 
approximate radius of gyration of a perfect sphere of ice at the 
simulation temperature. This parameter tends to zero for per- 
fect spheres. This measure again confirms that the crystalline 
clusters are reasonably spherical. 

We remark here that while these parameters measure spher- 
ical symmetry, they do not principally measure compactness: 
a cluster with small strings growing out of it on all sides will 
appear to be more spherically symmetric than a compact oval- 
shaped cluster, even though we might expect the oval cluster 
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to be more favourable (that is, having a lower surface area to 
bulk volume ratio). Nevertheless, we have not encountered 
any especially pathological clusters in our visual inspection of 
configurations. 



B. Nucleation pathway 

In order to compare further the nucleation behaviour we ob- 
serve to that of Quigley and Rodger, we have evaluated the 
global order parameters that they used along the nucleation 
pathways computed in this work. We used the Steinhardt- style 
26 an d Chau-Hard wick- style ' tetrahedrality parameters as 
defined by Quigley and Rodger, including their smoothing 
function, to ensure comparability of results. We thus have 



(a) 
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(8) 

where N is the number of particles, /(/*//) is the smoothing 
function, and we again use real spherical harmonics (the re- 
sults with complex spherical harmonics are, of course, identi- 
cal), and 87 

£ = ^ 1 1 1 Km) fir*) ifij ■ r ik + 1/3)) 2 , (9) 



i=\ j=\k>j 



where fij is the unit vector from particle i to particle j. The 
smoothing function is defined as 



f 1 



f(r) 



, cos 



(r/A-3.1)7l 
0.4 



1/2 







if r<3.1A, 

if 3.1A<r<3.5A, 

otherwise. 

(10) 

The bulk £ parameter for supercooled MW liquid water at 
220 K tends to ~0.21. For equilibrated cubic and hexagonal 
ice at 220 K, £ approaches 0.026, while 26 approaches 0.44 
for hexagonal and 0.51 for cubic ice. 

It is also important to note that as these global order pa- 
rameters measure the average order of the entire system, an 
ice cluster of a given size will produce different values of 26 
and £ depending on the overall size of the system. There- 
fore, as we wish to compare our results to Fig. 3 of Ref. 42, 
where the number of particles in the simulation was 576, we 
have calculated the global order parameters in such a way that 
we have only taken into account the nearest 576 particles from 
the centre of mass of the crystal nucleus as identified by our 
local order parameter. The mean values of 26 and £ along our 
pathway for the nucleation from a cubic ice seed are depicted 
in Fig. 6(a). At the top of our free energy barrier, we find 
26 ~ 0-10 and £ « 0.16. By contrast, Quigley and Rodger re- 
port that the saddle point of their two-dimensional free energy 
landscape occurs at 26 ~ 0-22 and £ « 0. 19. 42 Furthermore, in 
our simulations, the order parameters £ and 26 vary roughly 
linearly as the ice cluster grows, as would be expected for the 
growth of a crystal nucleus into a relatively unperturbed liq- 
uid. By contrast, Quigley and Rodger initially observe an in- 
crease in 26 with little change in £. The nucleation pathway 
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FIG. 6. (a) The global order parameters <26 an d C calculated along the 
local order parameter used to drive nucleation in our system. Error 
bars show the standard deviation for the population of configurations 
at each cluster size. The cyan solid line indicates the approximate 
location of the nucleation barrier as calculated by local order param- 
eters. The results depicted here refer to the 576 particles around the 
centre of mass of the ice nucleus for comparison with Ref. 42, and 
not to the full 1400-particle system, (b) An analogous diagram for a 
single brute-force trajectory at 21 OK. 



we observe therefore appears to differ significantly from the 
one reported by Quigley and Rodger. The sigmoidal form of 
their 26 -£ free energy landscape suggests their pathway in- 
volves an initial increase in the 26 orientational order of the 
whole system, not just of a crystal nucleus. This 'pre-ordering' 
then appears to make crystallisation more tractable on compu- 
tational timescales, even though this is not the natural pathway 
of nucleation. 

We have not been able to nucleate ice in a brute-force sim- 
ulation at 220 K; however, we assume that the nucleation be- 
haviour does not change significantly between this tempera- 
ture and other temperatures above the liquid transformation 
temperature of 202 K. An analogous diagram to Fig. 6(a) 
is shown for a representative brute-force simulation of nucle- 
ation at 210 K in Fig. 6(b). Even though the critical nucleus 
size is different, and the data are significantly more noisy, the 
basic linear behaviour of 26 an d C remains unchanged. We be- 
lieve that this brute-force behaviour is a good indication that 
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applying umbrella sampling does not change the natural nu- 
cleation pathway. 



C. Nucleation rate 

Configurations from the top of the barrier obtained in Monte 
Carlo umbrella sampling simulations were taken as starting 
points in isobaric-isothermal MD simulations (with a Nose- 
Hoover barostat and thermostat with relaxation times 2ps and 
lps, respectively, and a time step of 5fs). Velocity compo- 
nents were assigned from a normal distribution with mean 
zero and scaled to correspond to the simulation temperature 
of 220 K. 

For the system nucleating from a seed hexagonal crystal, 
we ran 40 simulations from distinct starting configurations at 
the top of the barrier. Of these, 17 exhibited nucleus growth 
and the remainder nucleus shrinkage. This roughly 1 : 1 ra- 
tio of growing and shrinking suggests that we are approxi- 
mately at the top of the barrier. From these simulations, we 
find the long-time attachment rate (averaged over all config- 
urations, not just the ones remaining close to the critical nu- 
cleus size 73 ) to be f Hcrit = 1.175 x 10 13 s _1 . The Zeldovich 
factor is Z 2 = 0.000344 and the number density of the liquid 
is Piiq = 3.3 x 10 28 m -3 . The probability of reaching the top 
of the free energy barrier is P(n c rit) = ex P [ — AG(n C n t )/fe^] = 
9.0 x 10~ 9 . The rate of MW ice nucleation at 220K is there- 
fore estimated to be R = 6.5 x 10 31 m -3 s _1 . 

We performed analogous simulations starting from the crit- 
ical cluster in the system which was grown directly from the 
supercooled liquid and has a critical nucleus size of 1 14. Of 
the 40 MD trajectories started from this point, 20 exhibited 
shrinking and 20 exhibited growth of the largest crystal nu- 
cleus. An equivalent analysis to that above gives a rate of 
nucleation of R = 2.6 x 10 29 m -3 s _1 . 

Using forward flux sampling with molecular 

dynamics, ' Li and co-workers have determined nu- 
cleation rates for the homogeneous nucleation of ice with the 
MW model ranging between 2.148 x 10 25 m" 3 s _1 at 220 K 
and 1.672 x 10~ 7 m~ 3 s -1 at 240 K, 5-10 orders of magnitude 
below the experimental values at the temperatures where these 
have been reported. Our results are somewhat higher than 
those reported by Li and co-workers. To an extent, we can 
rationalise this difference by noting that their use of forward 
flux sampling ensures a purely kinetically driven process; by 
contrast, our systems are locally equilibrated. This is not usu- 
ally a problem with algorithms such as forward flux sampling 
because local equilibration is generally considerably faster 
than progression to a higher value of the order parameter 
serving as a reaction co-ordinate; however, in systems such 
as water, local equilibration can be very slow. This serves to 
rationalise the difference in the calculated rate between the 
system we have grown directly from the supercooled liquid 
and the one simulated by Li and co-workers. 

A second issue is that there is no straightforward way for 
a simulation method to choose the pathway of fastest growth 
if this difference arises only for somewhat larger crystal nu- 
clei; indeed, the free energy profiles we have calculated look 
the same for small clusters. In forward flux sampling, the or- 



der parameter space is staged; however, if all the pathways in 
the initial stages are equally fast, then the pathways chosen in 
forward flux sampling will be the ones that are the dominant 
ones in numbers, not the ones that may, ultimately, prove to 
be faster. Once an ultimately slower pathway is chosen, the 
simulation is extremely unlikely to be able to backtrack and 
choose a faster pathway. We suggest that this is why Li and co- 
workers did not see the nucleation of a pure phase of ice, even 
though our results suggest that the rate of nucleation of pure 
hexagonal ice is two orders of magnitude higher than that of 
mixed ice. A greedy simulation algorithm such as forward flux 
sampling, despite sampling natural dynamics, is even more 
likely to proceed down the 'wrong' pathway than simulation 
methods which allow for local equilibration, such as umbrella 
sampling. Nevertheless, umbrella sampling is by no means 
immune to this problem: indeed, the hexagonal crystal that we 
have grown was only possible because we effectively imposed 
that nucleation pathway by seeding the system. 



IV. CONCLUSIONS 

We have simulated the homogeneous nucleation of ice as 
a function of the size of the largest crystalline cluster using 
the MW potential. The free energy barriers we have calcu- 
lated agree reasonably well with classical nucleation theory, 
and the interfacial free energies extracted from fits to the com- 
puted free energy profiles appear to be sensible. We have ob- 
served that the crystalline clusters grown directly from the su- 
percooled liquid tend to be cubic, although they have a higher 
free energy than those clusters grown from ideal crystalline 
seeds. 

A comparison with the free energy landscapes computed in 
the previous studies of homogeneous ice nucleation 36,42,43 at 
a similar undercooling (~20%), albeit with the TIP4P poten- 
tial, reveals that the free energy barrier calculated in this work 
is considerably lower. We suggest that previous simulations 
were not at equilibrium in path space; that is, the lowest free 
energy paths between the liquid and the crystal were not lo- 
cated. The high free energy barriers previously reported may 
be a result of non-local ' and non-rotationally-invariant 
order parameters used to track the progress of nucleation. We 
therefore hypothesise more generally that care should be ex- 
ercised when using global orientational order parameters, be- 
cause there may be a tendency for the pathways that are easiest 
to locate with rare event techniques not actually to be those 
with the lowest free energy, particularly for systems such as 
water, where crystal growth is slow. 

The simulations we have performed allow us to learn some 
lessons for future simulations of all-atom models of water. If 
the MW potential is a reasonable model for the energetics of 
water, as is strongly suggested by the work of Molinero and 
co-workers, 57,61 the difficulty in simulating water nucleation 
in all-atom models such as TIP4P may be more due to the 
slow kinetics of the process than an insurmountable free en- 
ergy barrier. We suggest that, to test this hypothesis, a study 
of the homogeneous nucleation of ice using an all-atom wa- 
ter potential and driven by a local order parameter would be a 
worthwhile, if difficult, endeavour. We are hopeful that the re- 
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suits for the thermodynamics of ice nucleation presented here 
will be a useful guide for such a study, which will in turn lead 
to a fuller understanding of the dynamics of ice nucleation. 
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